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Abstract. 

Context. Sticking of colliding dust particles through van der Waals forces is the first stage in the grain growth process in proto- 
planetary disks, eventually leading to the formation of comets, asteroids and planets. A key aspect of the collisional evolution 
is the coupling between dust and gas motions, which depends on the internal structure (porosity) of aggregates. 
Aims. To quantify the importance of the internal structure on the collisional evolution of particles, and to create a new coagu- 
lation model to investigate the difference between porous and compact coagulation in the context of a turbulent protoplanetary 
disk. 

Methods. We have developed simple prescriptions for the collisional evolution of porosity of grain-aggregates in grain-grain 
collisions. Three regimes can then be distinguished: 'hit-and-stick' at low velocities, with an increase in porosity; compaction 
at intermediate velocities, with a decrease of porosity; and fragmentation at high velocities. This study has been restricted to 
physical regimes where fragmentation is unimportant. The temporal evolution has been followed using a Monte Carlo coagula- 
tion code. 

Results. This collision model is applied to the conditions of the (gas dominated) protoplanetary disk, with an a T parameter 
characterising the turbulent viscosity. We can discern three different stages in the particle growth process. Initially, growth is 
driven by Brownian motion and the relatively low velocities involved lead to a rapid increase in porosity of the growing aggre- 
gate. The subsequent second stage is characterised by much higher, turbulent driven velocities and the particles compact. As 
they compact, their mass-to-surface area increases and eventually they enter the third stage, the settling out to the mid-plane. 
We find that when compared to standard, compact models of coagulation, porous growth delays the onset of settling, because 
the surface area-to-mass ratio is higher, a consequence of the build-up of porosity during the initial stages. As a result, particles 
grow orders of magnitudes larger in mass before they rain-out to the mid-plane. Depending on the precise value of ar T and on the 
position in the nebula, aggregates can grow to (porous) sizes of ~ 10 cm in a few thousand years. We also find that collisional 
energies are higher than in the limited PCA/CCA fractal models, thereby allowing aggregates to restructure. It is concluded that 
the microphysics of collisions plays a key role in the growth process. 

Key words. ISM: Dust - Planetary systems: formation, protoplanetary disks - Accretion disks 



1. Introduction 

Understanding the formation of planetary systems is one of 
the central themes of modern astrophysics. New stars form 
in molecular cloud cores when these cores contract under the 
influence of gravity. This contraction leads to the formatio n 
of a central object surrounded by a disk dShuetalJfl987l) . 
Planets are generally thought to form in these disks, but nei- 
ther the precise physical conditions required for, nor the pro- 
cesses involved in planetary body assemblage are well under- 
stood. Generally, it is thought that grain growth starts with 
the sticki ng of sub-micron-sized grains colliding at low ve- 
locities JWeidenschilling & Cuzzilll993l) . The sticking is then 
driven by weak van der Waals interaction forces between the 
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grains. Relative velocities may reflect Brownian motion or dif- 
ferences in coupling to turbulence in the disk. Eventually, when 
the grains have grown to ~cm-sizes, they will settle to the 
mid-plane of the disk, forming a thin dust lay er where f urther 
growth to planetesimal sizes c an take place JSafronovi ri 969 : 
Weidenschilli ng & Cuzzill993l) . 

There is much observational support for the growth of dust 
grains in protoplanetary disks from sub-micron to millime- 
tre size scales. In particular, the contrast of the 10 fim sili- 
cate emission feature relative to the local continuum shows 
that the grain size in the disk's photosphere - where these 
features originate - has increased from sub-micron sizes 
character istic for interstellar dust, to the m i cron-sized r ange 
Ivan Boekel etail l2003t iMeeus et alJ Eooi IPrzvgodda etatl 
2003e iKessler-Silacci et al.1 120061) . Furthermore, observations 
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of the continuum (sub)millimetre emission - originating from 
the deeper layers of protoplanetary disks - show typical 
grain sizes in the range of millimetres, outside the range 
of in terstellar grain sizes by many orders of magnitudes 
dBeckwith & Sargent! 199 lb . Additional support for the impor- 
tance of collisional grain growth follows from analytical stud- 
ies of solar system dust. In particular, many interplanetary dust 
particles (IDPs) collected at stratospheric altitudes, consist of a 
large number of small grains assemble d in a very open structure 
as expected for collisional aggregates dBrownleel 1 979h . These 
types of IDPs are thought to derive from comets and, indeed, 
comets may consist largely of such loosely bound aggregates. 
In addition, chondrules recovered from meteorites show dust 
rims which are generally attributed to collisional accretion pro- 
cesses in the solar neb ula before the meteorite and i ts parent 
body were assembled dMetzler et alll992tlCuzzil2004T) . 

The properties of the dust are of key importance to the evo- 
lution of protoplanetary disks. First and foremost, planet for- 
mation starts at the dust size and the dust characteristics at 
the smallest sizes will set the table for further growth. Second, 
the opacity is dominated by absorption and scattering by dust 
grains. Hence, the radiative transfer, temperature structure, as 
well as emission spectrum of protopla netary disks are con- 
troll ed by the characte ristics of the dust dBeckwith et all2 000: 
iBouwman et alJEo OO). Third, in turn, the temperature struc- 
ture will dominate the structure of the disk, including such as- 
pects as flaring. Fourth, dust grains provide convenient surfaces 
that can promote chemical reactions. Specifically, ice man- 
tles formed by accretion and reactions between simple pre- 
cursor molecules are wi despread in regions of star formation 
( Gib b et al1l2004l |Boogert et al.l Eo04). In fact, grain surfaces 
may be catalytically active in the warm gas of the inner d isks 
aroun d protostars, converting CO into CH4 dKress & Tielensl 
l200ll) . 

Most early studies of the coagulation process and the 
characteristics of the resulting aggregates assumed hit-and- 
stick collisions where randomly colliding partn ers stick at the 
point of initial c ontact jMeakinlfl985 IMeakin & Donni ri988: 
IOssenkop|l993l) . The structure of the aggregate then depends 
on whether the collision is between a cluster and a monomer 
(PCA) or between two clusters (CCA). The latter leads to very 
open and fluffy structures with fractal dimensions less than 2, 
while the former leads to more compact structures and a fr actal 
dimension (for large aggregates) near 3. Ossenkopl dl993l) also 
investigated the pre-fractal limit in which aggregates consist 
of < 1000 monomers. He provides simple analytical expres- 
sions for, e.g., the geometrical and collisional cross-section in 
the case of PCA and CCA aggregation. These expressions in- 
clude a structural parameter, x, which describes the openness 
(or fluffiness) of the particle. In this study we also introduce a 
structural parameter and confirm its importance in coagulation 
studies. 

Over the last decade much insight has been gained 
in the structure of collisional aggregat es through exten- 
sive, elegant, experimental studies (Blum - eUdJj2002l iBlum 
2004) supported by theoretical analysis dChokshiet aDl 19931 
Dominik & Tielens 19971) . These studies have revealed the im- 
portance of rolling of the constituent monomers for the result- 



ing aggregate structure. At low collision velocities, the hit-and- 
stick assumption is well justified but at high collision velocities, 
aggregates will absorb much of the collision energy by restruc- 
turing to a more compact configuration. At very high velocities, 
collisions will lead to disruption, fragmentation, and a decrease 
in particle mass. The critical velocities separating these colli- 
sional regimes are related to material properties such as surface 
free energy and Young's modulus as well as monomer size and 
the size of the clusters. 

While the porous and open structure of collisional aggre- 
gates is well recognised, their importance for the evolution of 
growing aggregates in a protoplanetary setting is largely ig- 
nored. Most theoretical studies repres ent growing aggregates 
either by an equivalent sphere (e.g., 1 Weidenschilling 1984a; 
Mizu no et al 1988; Tanaka et alJl2005t iNomura & Nakag awal 



2006) or adopt the fractal dimension linking mass and size 
chara c teristic for CCA or PCA growth (e.g., [W eidenschilling 



1997; 


Suttner & Yorke 


2001; 


Dullemond & Dominik 2005). 


Ossenkopfl d 19931) andlKempf et alJ 


( 1999) explicitly account 



for aggregates' internal structure in their numerical models, al- 
though, due to computational reasons, only a limited growth 
can be simulated. Indeed, the internal structure of collisional 
aggregates is the key to their subsequent growth. The coupling 
of aggregates to the turbulent motion of the gas is controlled 
by their surface area-to-mass ratio, while the relative velocity 
between the collision partners dictates in turn the restructuring 
of the resulting aggregate. Moreover, as a result of the growth 
process from sub-micron-sized monomers to cm-sized aggre- 
gates, the coupling to the gas velocity field may well change 
due to compaction. Indeed, compaction can be an important 
catalyst for aggregates to settle out in a mid-plane dust layer. 
Despite its importance for the collisional growth of aggregates 
in a protoplanetary environment, the evolution of porosity has 
not yet been theoretically investigated. The present paper fo- 
cuses precisely on this aspect of grain growth in protoplanetary 
disks. 

This paper is organised as follows. In Sect. |2]a model is 
presented for the treatment of the porosity as a dynamic vari- 
able. This entails defining how porosity, or rather the openness 
of aggregates, is related to the surface area-to-mass ratio (Sect. 
12. 21 as well as quantifying how it is affected by collisions (Sect. 
12. 31 . In Sect. |3 a Monte Carlo model is presented to compute 
the collisional evolution, taking full account of the collisional 
aspect and all features of the porosity model. Section @] then 
applies the model to the upper regions of the protoplanetary 
disks. Results from the porous model of this paper are com- 
pared to traditional, compact models. In Sect.|5]we discuss the 
differences of the new collision model with respect to PCA and 
CCA models of aggregation and also discuss our results from 
an observational perspective, before summarising the main re- 
sults in Sect. [6] 

2. Collision model 

Dust grains are dynamically coupled to the turbulent motions 
of the gas and 'suspended' in the (slowly accreting) proto- 
planetary nebula. Upon collisions, these small dust grains can 
stick due to van der Waals forces, forming larger aggregates. 
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Eventually, when these agglomerates have grown very large 
(~cm-sized), they will decouple from the gas motions and set- 
tle in a thin disk at the mid-plane. Further collisional growth 
in the mid-plane can then lead to the formation of planetesi- 
mals (~km-sized). Upto this point, growth is driven by weak 
van der Waals forces, but for km-sized planetesimals gravita- 
tional forces take over and rapid growth to a planet takes place. 
Here we focus on this process of small grains suspended in the 
nebula forming larger conglomerates. Coagulation is driven by 
the relative grain velocities. Velocities and the kinematics of 
dust in a turbulent nebula are discussed in Sect. 12.11 The fric- 
tional coupling between dust and gas depends largely on the 
area-to-mass ratio of the grains and hence on the internal struc- 
ture of the dust. Section l2~2l describes the relation between the 
area-to-mass ratio and the porosity of the dust agglomerates. 
In Sect. 12.31 we discuss experimental and theoretical studies 
on the microphysics of dust coagulation and develop a simple 
model, given in the form of easily applicable recipes, which 
describes how the mass and porosity of the dust evolve in col- 
lisions between two dust agglomerates. In Sect. 12.41 finally, 
we compare these recipes in the fractal limit to the frequently 
used formulations of Particle-Cluster Aggregation (PCA) and 
Cluster-Cluster Aggregation (CCA). 

2.1. The turbulent protoplanetary disk 

For the characterisation of the gas parameters of the proto- 
planetary disk we use th e minimum-mas s sol ar nebula (M SN) 
model as described by lHavashil Jl98ll) and iNakagawa et alJ 
Jl986l) . The surface gas density of the disk, S g , is assumed to 
fall off as a -1.5 power-law with heliocentric radius (R) and 
the temperature scales as R~ 1 ^ 2 . The vertical structure of the 
disk is assumed isothermal, resulting in a density distribution 
which is Gaussian in the height above the mid-plane, z. The 
scaleheight of the disk, H g , is an increasing function of radius, 
H g = Cg/Q. oc R 5 / 4 , with c g the local isothermal sound speed 
and Q the Keplerian rotation frequency. The thermal gas mo- 
tions will induce relative (Brownian) velocities between dust 
particles with a mean rms-relative velocity of 



Av Br ownian('«l,'W2) 



%kftT(mi + 1112) 
■Km\mi 



(1) 



with m\,m,2 the masses of the particles and Boltzmann's 
constant. Except for low mass particles of size < 10 yum, 
these velocities are negligibly small when compared to the 
relative velocities induced by the coupling with the turbulent 
gas. We assume that the turbulence is characterised by the 
Sha kura & Sunvaevl(ll973h ax-parameter, 



a T c z g /n 



(2) 



with Vq,€q an d h), respectively, the velocity, the size and the 
turn-over time of the largest eddies. According to the stan- 
dard (Kolmogorov) turbulence theory, turbulent energy is in- 
jected on the largest scales and transported to and eventually 
dissipated at the smallest eddies - characterised by turn-over 
time (or dissipation timescale), t s , velocity, v s , and scale size, 



f s , given by v m = v s € s , with v m the molecular viscosity. 1 We 
can then define the Reynolds number as, Re = vj/v m , giv- 
ing v s = Re~ 1 ^ 4 \o and f s = Re ~ l ^ 2 tn. If fa is assumed to b e 
(of the order of) the Kepler time dDubrulle & ValdettaroTT9 92). 
fa = 27r/Q~ I , the turbulence is fully characterised by the ctj- 
parameter (see Fig.[Q: 



fa = InQT 



1/2 

v = a T c g , 



e = a' 2 K 



t s = Re- 1/2 t 



v s =Re- l/ \ 



4 = Re~ 3f %. 



(3a) 



(3b) 



(3c) 



This specification of turbulence is of importance, for, to- 
gether with the friction times of two particles, it determines 
the (average root-mean-square) velocity between the particles, 
Avy, which in turn plays a crucial role in both the collision 
model of this section as well as in the time-evolution model of 
Sect. 12. 31 These relative velocities are a function of the friction 
time (Tf) of the particles - the time it takes a particle to react 
to changes in the motion of the surrounding gas - which in the 
Epstein limit is 2 



4cgp g A ' 



(4) 



where p g is the local mass density of the gas, m the mass of 
the particle and A its cross-section. In particular, if the friction 
time of a particle is much smaller than the turnover time of the 
smallest eddy (Tf <K f s ), the particle is coupled to all eddies and 
flows with the gas. Therefore, grains with Tf «c t s have highly 
correlated velocities. Eventually, however, due to growth and 
compaction, agglomerates will cross the Kolmogorov 'barrier' 
(i.e., Tf > f s ) and are then insensitive to eddies with turnover 
times smaller than Tf, since these eddies will have disappeared 
before they are capable of 'aligning' the particle's motion. At 
this point, grains can develop large relative motions (Fig.^. 

To calculate relative velocities accurately, con tributions 
from all eddies have to be included. IVoelk etal J lT98m started 
quantifying these effects by dividing eddies into several classes 
and subsequently integrated over the full Kolmogorov power 
spectrum. With the e xception of some specia l cases of the par- 
ticle's friction times dCuzzi & Hogadl2003l) in which Av can 
be expressed analytically, Av between pa rticles with differ- 
ent Tf can be presented only numerically JVoelk etal] [i980: 
Ifvlarkiewicz et alJll99lh . Iweidenschillingf ill 984 J hi) , however, 
has presented analytical fits, which have been frequentl y ap- 
plied in subsequent coagulation models (e.g.,|Su ttner & Yorkd 
200 ll iDullemond & DominiklEool iTanaka et alJl2005t) . For 



1 Vm = c eH m H/2p E cr m oi with jj and <x rao i, respectively, the mean 
molecular weight and mean molecular cross-section of the gas 
molecules. 

2 The Epstein limit holds for particles with sizes smaller than the 
mean-free-path of the gas, a < |/l ra f p . If this limit is exceeded, friction 
times increase by a factor ^fl//t m f D and quadratically s cale with radi us 
<Whipplell972t lmidenschil ling|l97^ISchraDler & Heriningfcooi) . 
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1/2 

Vq = Q' T ' C s 
[10 s cm s _1 ] 




• Ti = T 2 

■Tt > r 2 



Re 



1/2 



[3 X 10 3 s] 



[3 x 10' s] ti 

t = 2%jn 



Fig. 1. Sketch of the turbulence induced relative velocities, 
Av, as function of friction time, t\, for equal friction times, 
T\ - T2, (solid line) and different friction times T2 <k T2 
(dashed l ine) according to the analytic fits of Eqs. J5all5bl > (after 
IWeid enschilling 1984a). to and vo are set by the global proper- 
ties of the disk, while the range over which turbulence is impor- 
tant is determined by the Reynolds number, Re. Values between 
brackets denote typical values for to — 1 yr, c g = 10 5 cm s _1 
and a T = 10~ 4 . 

particle friction times, T\ and t 2 (j\ > t 2 ), less than to, these 
read, 3 



AV12 (Ti,T 2 ) 



(ji ~ t 2 ) 

i 

3 

! 1 + T 2 /T] "V t s 



n,r 2 < t s 



t s < T\ < to, 



(5a) 



where t\ is the larger of the two friction times. If t\ exceeds to 
the relative velocities become, 



Av 12 = • 



v 



vo- 



( Tl + T 2 )f() 

2t\t 2 



T 2 < t Q < Ti 



T\,T 2 > t . 



(5b) 



For example, in the regime where both friction times are small 
(ti , T2 < f s ) the turbulence induced relative velocity is negligi- 
ble when T2 ~ T\, but scales linearly with t\ when particle 2's 
mass-to-area ratio is much larger than that of particle 1 (Fig. 
0. Thus, in the t\,t 2 < t s regime particles with very differ- 
ent A/ra-ratios will preferentially collide, because differential 
velocities are then highest. When one of the particles' friction 
time exceeds t s the dependence on the absolute difference in Tf 
vanishes and the relative velocities now scale with the square 
root of m/A of the largest Tf. As can be seen in Fig.^relative 
velocities are rather insensitive to the precise ratio of the parti- 
cles' friction times in this regime. (Because of the simplicity of 
the expressions for Av in Eqs. ( Ball5bt the lines do not connect 



3 iDullemond & Dominild 120051) note that the second expression in 
Eq. i5al may not exceed v . 



at Ti = t s and T\ = to.) In the Tf > to regime, (Eq.|5bJ relative 
velocities would stop increasing and eventually become only a 
minor perturbation to the motion of the particle. These large r t 
regimes are, however, not reached in the early stages of coagu- 
lation considered in this paper. Here, it is clear that the relative 
grain velocities are regulated by the area-to-mass ratio of the 
colliding grains which sets the friction timescale (Eq. |4} and 
hence the velocities (Eas. l5al|5bl . 

2.2. Porosity of agglomerates 

For compact, solid particles, 4 the area-to-mass ratio scales lin- 
early with the size. However, for porous aggregates the Afm 
ratio depends on the internal structure of the aggregates. In this 
section, we discuss how the internal structure of the aggregates 
affects collisions and, conversely, how collisions affect the ag- 
gregates' internal structure. The internal structure is modelled 
using the enlargement parameter. Here, we define an enlarge- 
ment parameter if/ that is the ratio between its extended volume, 
V, and its compact volume, V*, i.e., 



V 

<A = — . 



lf/> 1. 



(6) 



V* is the volume the material occupies in its compacted state, 
i.e., when particles are packed most efficiently, and V the total 
(extended) volume of the particle. While V* reflects the mass of 
the particle, V is related to the spatial extent of the aggregate. 
Here, V has been defined as the volume corresponding to the 
geometric cross-section of the aggregate (see Fig. [3J». Figure|21 
shows three aggregates of 1 000 monomers, such that V* is the 
same, though with different degrees of compaction. The (from 
left to right) decreasing compaction translates to an increased 
geometric cross-section, A (outer circle), of the aggregates and 
hence to an increased iff. Using the linearity between V* and m 
and Eq. (|6j, two parameters, e.g., m and iff, determine A and, 
consequently, also determine the coupling with the gas (Eq.|4j. 

if/ — 1 then corresponds to the enlargement factor of 
a compact particle with specific density p s = m/V*. This 
does not necessarily correspond to a homogeneous particle (of 
zero porosity); it corresponds, however, to the lowest porosity 
that can be achieved by re-arranging the constituent particles 
(monomers) within the aggregate, without destroying this sub- 
structure through, e.g., melting. Since if/ > 1 we refer to if/ as 
the enlargement parameter or enlargement factor, if/ is related 
to the structural parameter, x, of lOssenkopfl \ 1 9931) . as x oc if/ 1 . 
Note, the close relationship between if/ and porosity; a higher 
if/ means more open aggregates, i.e., higher porosity. Therefore 
we will also use this more familiar reference for if/, but only in 
a qualitative sense. 

Since our enlargement parameter, if/, is defined in order to 
match the geometrical cross-section, A, of a particle, we refer 
to the resulting radii, a(if/), as geometrical. The other cross- 
section of importance in coagulation m odels is <x, the col- 
lisional cross-section. lOssenkopfl \ 19931) has consequentially 



4 The reader should note that the words 'particle', 'agglomerate' 
and 'aggregate' are frequently interchanged throughout this and other 
paragraphs. 
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porous size (V) 
compact size (V ) 




^6 



Fig. 2. Projections of fractal aggregates illustrating the relation between A, V and V*. All fractals consist of 1 000 monomers 
and have consequently the same compact volume, V* (inner circle). The black, outer circle gives the area, A, equal to the total 
projected surface area of the agglomerate. This area subsequently defines the volume, V (e.g., V ~ A 3 ^ 2 ) and the enlargement 
factor iff of the aggregate (Eq. |6j». (left) Compact aggregate, defining the compact volume V*. (centre) Porous aggregate. The 
geometric cross-section, A (outer circle), is larger than its compact equivalent (inner circle), (right) Even more porous aggregate. 
Arrows denote the compact and porous radii, a* and a. 



defined a 'toothing radius', a to oth, such that cr = 7r(fl t0 oth, l + 
fltooth, i) 2 ar, d provides expressions for a to oth f° r PCA and CCA 
aggregat es, cr is l arger t han A by a factor of order unity (see also 
iRrause & BlurrJ J2004)') and also depends slightly on the struc- 
ture of the aggregates, i.e., whether PCA or CCA. In this pa- 
per, however, we are mainly concerned with obtaining a model 
for iff and have therefore simply used a for the calculation of 
the collisional cross-section. Likewise, we have also ignored 
augmentations of cr due to e ffects such as, e.g ., charges and 
grain asy mmetries (see again lOssenkopJ fl 993)) and rotation 
of grains ( Pas zun & Dominikl200cj|) . 

The aggregates' internal structure has important conse- 
quences for the coagulation rate. The geometric cross-section 
scales, for example, as iff 2 / 3 . More subtle are the effects of iff on 
relative velocities, which, as seen above, depend critically on 
the A/ra-ratio of both aggregates. In coagulation models where 
grains are represented as compact bodies Tf increases linearly 
with size (A oc m 2 / 3 ;r t oc m 1 ^). However, in the initial stages 
of coagulation aggregates stick where they meet - a process 
characterised by a build-up of porous, fluffy structu res, which 
can be described by fractals. lMeakin & DonnN 19881) have com- 
puted the A/m ratio of an initially monodisperse population and 
found a profound flattening of this ratio as compared to com- 
pact models in which it decreases as m~'/ 3 . Often, in the 'hit- 
and-stick' regime the growth shows fractal behaviour and the 
cross-section can be directly parameterised as a power-law, i.e., 



Aocm 6 , §<<5<1. (7) 

The lower limit 8 = § corresponds to the growth of compact 
particles, whereas the upper limit of 5 = 1 denotes the aggrega- 
tion of chains or planar structures. In the cluster-cluster coag- 
ulation (CCA) model of lOssenkopil i ll 9931) 6 = Sc.ca = 0.95, a 
result that agrees well with the findings o flKempf et alJ (H"999). 
using N-particle simulations in the Brownian motion regime. 
More recent models, which als o include rotation of aggre- 
gates dPaszun & Doimnikll2006l) also confirm this exponent. 
The scaling of friction time and enlargement factor with mass 
then become 

Tfocm 1 " 5 , tffocmi 6 ' 1 , (8) 

such that for compact aggregation models (6 = 2/3) iff stays 
constant. On the other hand, in models of cluster-cluster aggre- 
gation (6 ~ 1), area scales roughly as mass and Tf stays constant 
or is only weakly dependent on mass, while iff increases mono- 
tonically. Thus, in CCA relative velocities are rather insensitive 
to growth. As a result, collisions are also less energetic in CCA 
models, e.g., as compared to compact aggregation models. 

The key to the coagulation process in protoplanetary disks 
is the coupling of the dust to the turbulent motions of the gas 
and the resulting velocity distribution. In essence, the enlarge- 
ment parameter iff provides a relationship between mass and 
surface area for growing aggregates which controls this gas- 
dust coupling. Equation provides a relation for the evolution 
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of if/, but this relation is only valid in certain specific aggre- 
gation models, e.g., CCA-coagulation, where similar, equally 
sized aggregates meet. In reality, however, collisions between 
particles of all kinds of sizes will occur, although, dependent on 
the parameters that determine Av, some are just more probable 
than others. In the end, the growth of grains in protoplanetary 
disks is controlled by the individual collisions between two ag- 
gregates. Therefore, we have to provide prescriptions for the 
outcome of all possible collisional encounters, i.e., all relevant 
combinations of m, iff and Av. 

2.3. The collision model 

We consider a collision between two particles with the aim of 
applying the results to a true coagulation model. Essentially, 
we have to provide a recipe for the enlargement factor, after 
the collision. Two relevant parameter sets enter into the new i/r: 
the progenitor masses and enlargement factors (i.e., the (ot,, iffi) 
of the colliding particles). Moreover, the collision energy, 



radius [cm] 



1 m x m 2 2 
E = Av 

2 mi + ni2 



(9) 



with /j. the reduced mass, is of key importance. At very low 
velocities, collisions between two aggregates will lead to stick- 
ing without internal restructuring, i.e., the particles will stick 
where they make first contact. At moderate velocities, the in- 
ternal structure of the resulting aggregate will adjust - dissi- 
pating a major fraction of the kinetic collision energy - re- 
sulting in a compaction of the aggregates. Finally, at very high 
collision velocities, the colliding aggregates will fragment into 
smaller units and this ca n lead to su bstantial erosion. Following 
the numerical model of Dominik & Tielens ( 1997) and the ex- 
perimental studies bv lBlum & Wurm (2000), these collisional 
regimes are distinguished by the following critical (collision) 
energies: 

- £restr = 5E ro n, the energy needed for the onset of com- 
paction; 

- £max-c — N c ■ E IO fi, the energy at which aggregates are 
maximum compressed. Here, N c is the total number of con- 
tact surfaces (between monomers) of the agglomerate. For 
a very open, fluffy agglomerate, N c = N. With increasing 
compaction the number of contacts will increase and for a 
cubic packing N c = 3N. Here, for simplicity, we will adopt 
N c = N. 

- £frag — 3A^ C • £)>reak, the energy needed for (the onset of) 
fragmentation of the agglomerate. Here £break is the energy 
required to break a bond between two monomers. Its mag- 
nitude is of similar order as E m \\. 

For monomers of the same size E m \\ is given by 
JPominik & Tielensll997llBlum & Wurml20 00) 



(10) 



with «o the radius of the monomer, y the specific surface 
adhesion energy and £ cl -i t the yield displacement at which 
rolli ng commence s. F ro u, the rolling force, was measured by 
iHeim et alJ £999) to be F mn = (8.5 ± 1.6) x 10~ 5 dyn for 
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Fig. 3. The collision regimes as function of total particle 
mass and relative velocity. Thick dashed lines indicate the 
critical energies for the onset of rolling and fragmentation. 
Parameters are that of quartz particles (p s = 3.0 g cm~ 3 ,y = 
25.0 ergs cm -2 ,ao = 0.1 fim) and for the Ef rdg line we have 
assumed equal masses (mi = m^) and Ebreak = £roii- Arrows 
indicate how the critical energy lines shift with increasing 
monomer size and mass-ratio, m\lm.2. 

uncoated SiOa-spheres of surface energy density y = 14 ± 
2 ergs crrT 2 . We adopt this value for F ro n and assume pro- 
portionality with y when applying it to other materials. The 
monomer size, ao, is also an important model parameter di- 
rectly affecting the outcome of the collision at a given mass; 
a smaller «o provides more contacts (for the same mass) and 
the agglomerate is more resistant to breakup. With these en- 
ergy thresholds, three qualitatively different collision regimes 
can then be discerned (see Fig.[3J: 

i E < E lestr : No internal restructuring. The particles stick 
where they meet ('hit-and-stick' growt h) as in t raditional, 
lattice-based, aggregation models (e.g.. iMeakinll 19881) - a 
process leading to fractal aggregates. 

ii Erestr < E < Ef rds : Restructuring (compaction) of aggre- 
gates. 

Hi E > Zifrag- Initially, loss of monomers, but for high ener- 
gies complete fragmentation (e.g., catastrophic collision). 
This phase requires a recipe for the mass distribution of the 
fragments. In this paper, fragmenting collisions are avoided 
by, e.g., 'choosing' moderate aj and stopping coagulation 
for particles that reach a critical friction time (see Sect.0}. 
Within these constraints, our results are therefore not com- 
promised by ignoring the fragmentation issue. 

From Fig. [3] it is clear that, when starting with small parti- 
cles, growth will initially be in the fractal regime. This frac- 
tal growth will be followed by a period in which the colli- 
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sions will promote the restructuring of the growing agglom- 
erate. Fragmentation becomes only important for velocities in 
excess of 10 3 cm s _1 . Since we assume that every contact can 
absorb a unit energy E ro u, the Ef mg line in Figgis independent 
of total mass. The Av a 10 3 cm s _1 limit translates to a critical 
value for the turbulent aj parameter: ctj < 1CT 2 . Within Fig. [3] 
the precise 'growth path' of the agglomerate will be controlled 
by evolution of the relative velocities and hence A jm or, equiv- 
alently, tfr. We will now construct recipes for iff in, respectively, 
i) the fractal and ii) the compaction regime. 



2.3.1. E < E n 



hit-and-stick 



Besides the usual CCA and PCA formalisms, there have 
been a few attempts to give prescriptions for the evolving 
inter nal structure of aggregates in th e hit-and-stick regime. 
iKostoglou & Kons tandopoulos (2001) discuss a formalism for 
obtaining the new fractal dimension in terms of the sizes and 
fractal dimensions of the two colliding aggregates. One point 
is, however, that, apart from the fractal dimension, another pa- 
rameter - the prefactor - is needed to fully describe the fr actal, 
although it is usually of order unity. lOssenkonfl il 19931) . like 
this study, introduces only one structural parameter and inter- 
polates between the CCA and PCA limits. We will also follow 
this idea, but use a different interpolation mechanism. 

We recognise that in the pure sticking regime most col- 
lisions are between evolved, fluffy aggregates, since the size 
distribution evolves progressively toward larger sizes. For low 
velocity-mass combinations (Fig. [3}, where restructuring is 
unimportant, the collisional growth then resembles the CCA 
growth process the most. We therefore simply rewrite the frac- 
tal law in terms of the individual masses of the particles and 
keep the CCA exponent, 



m\ ) 



mi > nil- 



(ID 



Although collisions between different particles are included 
in Eq. dl It . we still adopt the CCA-characteristic exponent 
(<5cca = 0.95) to ensure that for the 'pure' CCA case (mi = 
ni2\ Ai = A2) this prescripti on is in accordance with d etailed 
numerical studies dMeakin & Donnl Il988l lOssenkopJ Il993l 
lKempfetalJll999t IPaszun & Dominikll2006l) . There is. how- 
ever, a modification to Eq. (lilt that must be made. The term 
in brackets in Eq. (II It determines the amount of increase in 
A in the fractal regime. Because fractal growth results from 
inefficient packing of voluminous objects, it is clear that this 
term must include parameters describing the spatial extent of 
the collision partners. These cannot, however, be given by the 
masses of the particles, since m (alone) does not reflect a spa- 
tial size. For example, if we would replace one of the aggre- 
gates by one of the same mass, but of lower porosity (i.e., a 
more compact aggregate), its volume is obviously smaller and 
packing becomes more efficient. These effects, however, are not 
reflected in Eq. (II It . For these reasons, we replace m by V in 
Eq. (HD, 



A=Ai 



V1+V2 
Vi 



#CCA 



Vi > v 2 . 



(12) 



Note that for particles of the same internal density (porosity) 
Eq. H2\ and Eq. (II Q agree, such that Eq. H2\ also can be seen 
as an extrapolation from the CCA case, but one that takes ac- 
count of the different internal structures of the collision part- 
ners. Using the relation A oc y 2 / 3 and Eq. (|6}, Eq. Jl 2i can be 
re-written in terms of m and iff only 



<A = Mm 1 + 



ni2i]/2 
muffi 



4 (V. 



(13) 



with {ijj) m the mass-averaged enlargement factor of the colli 
sion partners, 

mnj/i + ni2if/2 



(4>)m = 



m\ + nil 



(14) 



In CCA coagulation (mi = mi and if/i = 1//2) we recover Eq. 

but Eq. Jl 3i now includes all collisions in the hit-and-stick 
regime. For example, if a large, fluffy aggregate sticks to a com- 
pact one, the enlargement factor of the newly formed aggre- 
gate is higher than the mass-averaged enlargement factor of the 
progenitor particles, (t^) m , but smaller than that of the fluffy 
collision partner. In Sect. 12.41 it will be shown, however, that 
the {ip) m term underestimates the porous growth when one of 
the particles is very small, i.e., in PCA-like collisions. This is 
solved by adding to Eq. (II 31 a term that compensates for these 
cases and our final recipe then becomes 



niiij/ij 



§<5cca-1 



(15) 



where i^add, a term important only for small m%, is explained in 
Sect. 1231 

2.3.2. E > -EVestr.: compaction 

In the compaction limit monomers are restructured at the ex- 
pense of the porou s volume of the aggregate s. Following the 
theoretical study of Dominik & Tielensl I I 19971) we will assume 
that the (relative) amount of compaction, AV p /V p , scales lin- 
early with collisional energy, i.e., AV p /V p = -fc = E/E max . c = 
—E/(N ■ E ro ii), where V p = V - V* denotes the porous volume 
within the aggregate and N is the total number of monomers 
present. 5 Essentially, this implies that the collision energy is 
used to set individual monomers in an agglomerate rolling and 
that this rolling is only stopped when an additional contact is 
made, resulting in compaction. Recalling that V = iffV* and 
V p = V - V* = (if/ - 1) V* with V* proportional to mass, we then 
find that the porous volume after colliding is 

(V* + V 2 *)(iA - 1) = (1 - fc) [ W, - 1) + V 2 *0A2 - D] ■ (16) 
And the new enlargement factor 



<A-1 = d-/c)- 



1 



mi + m>2 

= (l-/ C )(Wm-l), 



•(tojOAi - l)+m 2 (^ 2 -l)) 



(17) 



Note that we start compaction already at E = E mtt instead of E = 
5£ I0 u. We have found, however, that the simulations are insensitive to 
the precise energy at which compaction starts. 
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Fig. 4. Relative compaction as function of the mass ratio, 
1112/1111, at various collisional velocities. The expression on the 
y-axis measures the compaction relative to particle 1 and is 
a function of mass ratio, m%lm\, only (see Eg. I17>. Here if/i 
is the enlargement factor of the most massive aggregate, i.e., 
m\ > ni2, and plots are shown for fa ^ <Ai (grey lines) and 
<A2 = 4>\ (black lines). At low mass ratios the curves converge. 
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Fig. 5. The enlargement factor in the PCA and CCA limits 
as a function of the number of mono mers (AO of the p arti- 
cle. The thick grey lines show the fits of Ossen kopfl dl993h for 
CCA (top) and PCA (bottom) coagulation. The solid line corre- 
sponds to the limiting cases of Eq. Jl 3i . The dashed lines show 
the same limits, but now the ze ro point lies at N — 40 (after 
which the PCA/CCA curves of lOssenkopJ dl993l) diverge in 
porosity) and includes the if/^d correction term (Eq.fT^i. 



with (if/) m again the mass-averaged enlargement factor of the 
(two) collision partners (Ea.ll4>. We illustrate Eq. (I17> in Fig. 
|4]for the limiting cases of equal porosity collisions (if/i = \f/2\ 
black lines) and very porous vs. compact particle collisions 
» \j/2\ grey lines). Higher mass-ratios give higher colli- 
sion energies (higher /j) and hence more compaction. If veloc- 
ities are low, (v < 100 cm s -1 ) the net compaction occurs pri- 
marily through the (if/) m factor and the curves converge on the 
cm s (thick) line. Only when v > 100 cm s _1 does the /ce- 
faclor start to become important. Collisions at velocities higher 
than 1 000 cm s can result in fragmentation if the mass-ratios 
are similar. For low mass ratios the (if/) m is determined by xf/y 
(the highest mass) and there is little difference between the two 
limiting cases. 

2.4. Porosity increase in the PCA and CCA limits 

In the discussion on the if/ recipes in Sect. l2.3.T1 we have used 
the CCA fractal exponent (6qca = 0.95) as the starting point 
and extrapolated this empirical relation to Eq. dl 3I >: a general 
recipe for all collisions in the 'hit-and-stick' regime. By defi- 
nition, CCA-collisions are incorporated into this recipe and we 
may expect Eq. Jl 3i also to account for collisions between ag- 
gregates having about the same size. The PCA model, where 
one monomer collides with an agglomerate, is the opposite ex- 
treme. If we take the PCA-limit of Eq. dl 3b . i.e., m\ » ;«2 ~ 
mo, with mo the monomer mass, the change in ip after addition 



of a monomer becomes 

ifr N +i ~i[>N + ^r (§<5cca<A2 -fa); N » N 2 , (18) 

with N - Ni - ml mo the number of monomers of the PCA 
agglomerate and Ni — 1 for monomers. Thus, \fij^ goes to 
§<5cca<A2 ~ 1.5 in the PCA asymptotic limit (N 2 = 1; ife = I)- 6 
The fact that an asymptotic limit is reached can be under- 
stood intuitively, since there must be a point at which the in- 
ward penetration of monomers into the centre of the aggre- 
gate, which decreases iff, starts to balance the porous growth 
due to hit-and-stick collisions at the surface. The asymptotic 
limit of iff x 1.5, however, is much lower than typical PCA 
models indicate (iff * 10) as is illustrated in Fig. [5] where the 
PCA/CCA limits of our model (solid li nes) are co mpared to de- 
tailed numerical simulations of Ossenkopf ( 1993) (thick lines). 
Equation Jl 3I > thus underpredicts the porous growth for PCA- 
like collisions in which one of the particles is small; a result 
not too surprising since it originates from the CCA fractal law 
(Eq. II 11 . which is constructed to apply only for similar (i.e., 
equal-sized) particles. For these reasons we add to Eq. Jl 3i a 
term that compensates the -Niifj^lN term in Eq. Jl 8i . 

(Aadd = B — if/i exp [-yu/m F ] , (19) 
mi 

6 Here we take tf/i = 1 as the enlargement factor of single 
monomers. 
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where the exponential ensures if/^ is unimportant for colli- 
sions between particles well above a certain mass-scale, m$. 
With B = 1.0 and m P _ 10 mn we f ind good correspondence 
with the results of Ossenkorjf (1993). In Fig.|5]the new CCA 
and PCA limits are shown by the dashed curves, where we have 
shifted the 'zero point' from N - 1 to N = 40, i.e, the start- 
ing point for if/ after which we recursively apply Eq. J15I . The 
transition toward fractal behaviour emerges only after this point 
and w e therefore directly take if/ from the results of Ossenkopf 
ill 9931) when N < 40. 

Although, with Eq. ( 1191 for if/ d dd, Eq. (I15> does achieve 
the right PCA/CCA fractal limits, we do not claim they actu- 
ally provide a model for if/. The collision recipes are based on 
empirical findings and extrapolations from these. However, in 
contrast to the CCA and PCA limiting collisional growth mod- 
els, where if/ can be directly parameterised in a single expo- 
nent, we recognise that the internal structure of the aggregates 
is changed during - and caused by - the collisional growth pro- 
cess. This is the main qualitative difference of our collision 
model captured in Eq. (1151 . This equation, together with Eq. 
d!7i . provides recipes for the collisional evolution of the en- 
largement factor, which can be easily incorporated into time- 
dependent coagulation models. We acknowledge this is an ac- 
tive area of research and future model efforts may well improve 
on the present formulation. 

3. Monte Carlo Coagulation 

3.1. Outline 

To determine the implications the new collision model has for 
the solar nebula, e.g., as compared to collision models where 
mass is the only parameter, it must be embedded in a coag- 
ulation model that evolves the particle distribution function, 
f(x). f(x, t) gives the number density of particles with a set 
of properties (parameters) {jc,} at time t. In compact coagula- 
tion models all properties depend only on mass, so fix, t) = 
f(m, t). In the model described in Sect. |3 however, the parti- 
cle's enlargement factor has been included as an independent 
parameter such that / becomes a function of three variables, 
/(x, f) = f(m, if/, t). The coagulation equation which describes 
the evolution of /(x) is 

df(x,t) 



dt 



= -/(*, 



j dx'f(x',t)K(x,x')+ 



(20) 



+ I j dx'dx"/(x', t)f(x", t)K(x',x", t) <J k (x - T (x', x"))) , 

with K = crAv the collision kernel, 6^ the Kronecker 5-function 
and T the collision function, which maps in the case of stick- 
ing 2n parameters (those of x' and x") to n, where n is the 
number of independent parameters with which a particle is 
characterised. Equation d20l o f course is just an ext ension of 
the Smoluchowski equation 7 ( Smolu chowskil \l9 1 6l) to more 



7 Which reads: 
df(m) 
dt 



f(m) 



J dm' K(m,m')f(m')+ 

If,, , 

— I dm K(m ,m — m )f(m )f(m - in), 



than one dimension. Applied to the formalism in Sect.[5]it de- 
scribes the loss of particles in state x = (m, if/) due to colli- 
sions with any other particle (first term) and the gain of 'x- 
particles' that happen to be formed out of any suitable col- 
lision between two other particles (second term). Applied to 
the findings in Sect. |3 T symbolises the collision recipes with 
F(mi ,if/\,m2, if/2) - (m\ + m-i, if/) and if/ is given by Eq. d!5l > or 
Eq. d!7i . dependent on the collisional energy. 

One approach to implement coagulation is to numeri- 
cally integrate the Smoluchowski equation. However, it is 
immediately clear that numerically integrating Eq. d20l be- 
comes a daunting exercise. Integrating the ordinary (1- 
dimensional) Smoluchowski equation is already a non-trivial 
matter. Problems of near cancellation (the gain terms often 
equal the loss terms), mass conservation (systematic propa- 
gation of errors) and the problem concer ning the determina- 
tion of a time-step must all be tackled. ( Dullemond & Dominik 
(2005) in their Appendix B give an overview of the subtleties 
involved.) The F-factor in Eq. d20i gives a further complication 
since there is no such thing as 'conservation of porosity' and 
the 5-factor cannot be easily integrated away. Although there is 
no fundam e ntal re ason against the binning method - see, e.g., 
lOssenkopH Jl993h who solves the Smoluchowski equation in 
two dimensions - these issues make the whole procedure quite 
elaborate. We felt that much of the simplicity of the collision 
model of Sect.|2]would be 'buried' by numerical integration of 
a 2d-Smoluchowski equation and therefore have found it suit- 
able to opt for an approach using direct Monte-Carlo simula- 
tion (DSMC) techniques. 

The simplicity of using Monte-Carlo methods for coag- 
ulation problems is appealing. Basically, TV-particles are dis- 
tributed over a volume *V (see Fig. |6j\). The evolution then 
boils down to the determination of the two particles which 
are involved in the next collision. We hereafter assume that 
the particles are well mixed, i.e., no potential is present, such 
that the determination of the next collision is governed by 
basic stochastic principles. Then the probability of a colli- 
sion between particles i and j is given by the collision rate, 
Cij - KjjfV in which K,j is the collision kernel. A random 
number determines which of the jN(N - 1) possible collisions 
will be the next. The collision then creates a new particle, af- 
ter which the {Cy} must be updated and the procedure repeats 
itself. 

The advantages of such an approach are obvious. Most 
striking perhaps is the 'physical character' of Monte-Carlo sim- 
ulations. The growth-evolution of individual particles is di- 
rectly traced and can be studied. The algorithm does not use 
the distribution function, /, in a direct way; it is obtained indi- 
rectly by binning the particles. Secondly, the above described 
method is exact, i.e., no 'time vs. accuracy' considerations have 
to be made in choosing the time-step At; instead, At - the inter- 
collision time - is an outcome of the stochastic coagulation pro- 



(21) 



describing losses of m due to all collisions with m (first term on right 
hand side) and gains in the distribution of m due to collisions between 
m' and m — m' (second ter m), where t he fac tor j ensures collisions are 
not twice accounted for. Ossenkoof 1 1993) provides a general exten- 
sion of the Smoluchowski equation including source and sink terms. 
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Fig. 6. Illustration of the adopted Monte Carlo technique. (A) Af -particles are assumed to be uniformly distributed in a box of 
volume < V. The collision rate between particle i and j is Cy = cryAvy/'V. A random number determines which two particles will 
collide. (B) After the collision the number of particles is restored by duplicating one of the remaining N—l particles. Physically, 
this can be interpreted as an expansion of *V (heavily exaggerated in this figure) to a volume at which *V contains Af -particles 
again. The hypothesis of this method is that the collisional evolution within *V is representative for the coagulation of the total 
space under consideration. 



cess as it is in nature. Furthermore, due to its stochastic nature, 
no Monte-Carlo simulation is the same. A series of (indepen- 
dent) runs gives at once a measure of the statistical spread in the 
distribution. Note that the fluctuations around the average are 
a combination of real stochastic noise and random noise, but 
it is qualitatively different from the Smoluchowski approach, 
which describes the evolution of the mean of all possibl e re- 
alisations and is therefore completely deterministic (Gillespie 
119771) . From a practical point of view the straightforwardness 
of the DSMC-method makes that there is no need for resorting 
to 'control parameters' like those required in the numerical- 
integration method. 

The DSMC-method, however, has its limitations. It can be 
immediately seen, for instance, that when having started with 
N (say identical) particles of mass mo in a fixed volume, these 
will over time pile up in one agglomerate of mfl na i = Nm^. 
The accuracy then steadily decreases during the simulation 
(in MC-simulations the statistical error scales proportional to 
N~ 1 ^ 2 ) and most of the computing time is spent in the first 
few (quite uninteresting) collisions. Consequently, to achieve 
orders of magnitude growth, the initial number of particles 
must be extremely large. And because the calculation of the 
{Cy} goes proportional to A^ 2 (every particle can collide with 
each other), it becomes clear that this method becomes im- 
practicable. To counter the dependence on large initial parti- 
cle numbers, one can also try to preserve the number of parti- 
cles during the simulation. This can be done, for instance, by 
'tossing-up' the particle-distributi on when the number of parti- 
cles reaches N/2 as described by Liffman ( 1991 ). A more nat - 
ural approach perhaps, given by Smith & Mat soukasl (|l998), 



and adopted here, keeps the number of particles constant at 
each step; every time a collision takes place one of the remain- 
ing N - I particles is randomly duplicated such that the num- 
ber of particles throughout the simulation stays the same. This 
procedure can graphically be represented as an increase in the 
simulated volume, 'V (see Fig.|6j}), under the assumption that 
the c ollisiona l evolution ou tside of *V progresses identically. 
Smith & Matsoukas (1998) have shown that the error in the 
coagulation-scheme now scales logarithmically with the extent 
of growth, - or growth factor (GF), defined here as the mean 
mass over the initial mean mass of the population - much im- 
proved over the constant-'V case, where the error has a square- 
root dependence on GF. We might worry though about the con- 
sequences of the duplication mechanism. It causes a certain de- 
gree of 'inbreeding', which effects we cannot quantify directly. 
ISmith & "M atsoukas ( 1998) show it is unimportant for the con- 
stant or Brownian kernels used in their studies. However, these 
kernels are known to behave gently, i.e., they are rather insen- 
sitive to irregular changes in the population since the varia- 
tions in the Ky are small. Perhaps, more erratic kernels are 
more sensitive to the 'duplication mechanism', but we might 
equally well attribute this sensitivity to the stochastic nature of 
the coagulation process. At any rate, these consequences are 
best quantified by running the code multiple times. 



3.2. Implementation 

In implementing the DSMC approach we follow the 'full con- 
ditioning method' of Gillespie! (119771) . This involves the calcu- 
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lation and updating of partial sums, C,, 

N 

C,= ^C, 7 , i = l,...,N-l, (22) 

j=i+l 

with Cy = Kjj/'V = o-jj(a i ,aj)Av i j(T i ,Tj)/ r V(K) (here k is the 
total number of collisions since the start of the simulation). 
These N - I quantities are stored in the memory of the com- 
puter. Got = Tji Cj is the total coagulation rate. The probability 
density function, P(t, i, j), i.e., the probability that the next col- 
lision will occur in time-interval (f , t +dt) and involve s particles 
i and j (/ < j) can then be written as ( Gilles pie! 19771) 

P(t,i,j) = Cyexp[-Ctotfl 

= (C t0 , exp [-C tot f] ) x (Cj/Cw) x (Cij/d). (23) 

Three random deviates, r, s [0,1], then determine succes- 
sively: T) the time it takes until the next collision takes place, 
t = C to | ln(l/ri); //) the first particle (z) to collide, by summing 
over the C,-s (starting with i = 1) until r2C tot is exceeded (this 
fixes /); and ///) its collision partner (J) by summing the Cy-s 
over the /-in dex (starting w ith j = i + 1) until the value r$Ci 
is exceeded (Gillespie 1977, Eq. 19). The outcome of the col- 
lision is evaluated using the relevant equations in Sect. |2 and 
the new particle is stored in the /-slot. Another random num- 
ber then determines which of the N - I particles (excluding j) 
will be duplicated and this one is stored in the j'-slot. Having 
created, removed and duplicated particles, all of the C,-s need 
to be updated. This implies only the subtraction/addition of the 
Cy-s that have changed, not the re-computation of Eq. i22\ . 
Moreover, the 'duplication procedure' entails a rescaling of the 
simulated volume, < V, such that the spatial density of solids, 
Pd = Yui mil'Vix) remains constant. The algorithm can then be 
repeated. All these steps are order- /V calculations at worst; most 
time-consuming are the determination of j (for which Cy has 
to be calculated) and the update of the {C,}. To achieve a given 
GF another factor /V in computation time is needed, 8 bringing 
the total CPU-time proportional to /V 2 . These procedures are 
graphically summarised in Fig.0 

It is possible, however, to save some CPU time by taking 
the duplicates together in the calculation of the collision rates. 
If there are (g, - 1) copies of particle /, it would be a waste of 
time to calculate the (same) rates g, times. Rather, gi can be 
absorbed in the calculation of the (combined) coagulation rate. 
Cy = gig jdj is then the rate at which one of the /-particles col- 
lides with one of the j'-particles (/ + j) and C„ = jgiigi - 1)C„ 
between duplicates. The CPU time per step is now proportional 
to N s , the total number of distinct particles, i.e., excluding du- 
plicates, with 2, = ! i gi = N. To think of it in biological terms: 
gi gives the multiplicity of species /; N s the total number of 
species; and /V the total number of living creatures. 

8 Due to the duplication, the mean mass of the system increases 
with a factor (N + l)/N. The growth factor after /(-steps then becomes 

(N+lY 

<* = (—). (24) 
Thus, InGF = /cln(l + A^ 1 ) ss k/N if N » 1 and k « AHnGF. 



Calculate/Adjust collision rates: 

Ky = AVjjO-jj : Collision kernel 
Cjj = gjgjKjj / V : Collision rates 
Cj; C tot : Partial, total rates 

Find collision pair (i,j) 

i = Func ( {C|}, r) 
j = Func ( Cj .{Cy}, r) 
At = - In ( Co.- r ) 



Perform collision 

Determine new (m, ip) 
Have particles rained-out? 

Pick duplicates 

Preserve either: 

N : cnst #particles 
OR N s : cnst CPU-time per cycle 
OR N x N s : cnst growth per cycle 

Update parameters: 

t = t + Al 
Add duplicates 

Determine new V [keep p cnst] 
Adjust collision rates 

I 

Fig. 7. Flow chart of the MC-coagulation method. One cycle 
corresponds to one collision. 'Func' like in / = Func( {C,}, r ) 
indicates that i is a function of the C, values and a random de- 
viate, r. The procedure is further explained in the text. 

3.3. Tests 

The Monte-Carlo coagulation model described above is tested 
with kernels that have an analytic solution. These are I) the 
constant kernel, K,j = 1, and //) the linear kernel, K,j = |(m,- + 
m.j). The evolutio n of the mean mass of the distribution, (m) for 
these kernels is JSilk & Takahashill979tlTanaka & Nakazawa 
11994 

|too(1 + hi) constant kernel 

(m) = \ r j -i (25) 

I mo exp \jt\ linear kernel, 

where the distribution at t = is monodisperse of mass 
mo- Well-known coagulation models have either K oc m 1 ^ 
(Brownian coagulation) or K oc m 2 ^ 3 (geometric area), but 
here, due to changes in if/ and Av, we should be prepared for 
K to vary with time. Thus, it is important that the Monte-Carlo 
code passes both these tests. Initial conditions for these test 
cases are monodisperse with all relevant parameters put at 1 

9 The mean mass of the population is inversely proportional to the 
number of particles per unit volume. 
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Fig. 8. Test of the Monte Carlo coagulation code - constant 
kernel (Kn = 1), 20000 particles. At (dimensionless) times 
t = 1, 10 2 , 10 4 , 10 6 , 10 8 particles are binned and the distribution 
function / is computed (symbols). The analytical solutions at 
these times are overplotted by the solid lines. The error bars 
(hardly visible) show the spread averaged over 10 runs. The 
dotted lines show the distribution function if all the bins would 
be occupied by only 1 particle - it is an auxiliary line with 
slope 1 showing the lower limit of the (individual) distribution 
function. 

at t = (i.e, m - 1 and f(t = 0) = 1, p<j = 1, p s = 1) 
and we do not take porosity into account (iff — 1 always). At 
various times the particles are binned by mass and the dis- 
tribution function f(m) is determined by summing over the 
masses in the bin and dividing by the width of the bin (to 
get the spectrum) and the volume of the simulation (to get 
the density). Multiple runs of the simulation then provide the 
spread in /. Figures [8] and |9]present the results. On the y-axis 
f(m) is multiplied by m 2 to show the mass-density per loga- 
rithmic bin. Analytical solutions (Ta naka & Nakazawal Tl994) 
are overplotted by solid curves, while the dotted line shows 
the (hypothetical) distribution function if all the bins would be 
occupied by 1 particle. Thus, the dotted line corresponds to 
m 2 f(m) = m 1 fVwb ~ mfV, because the widths of the bins, 
w>b, are also exponentially distributed. In a single simulation, 
the distribution function should lie above this (auxiliary) line 
and the vertical distance to this line is a measure for the num- 
ber of particles in a bin. 

Figure|S]shows that the code passes the constant kernel test 
with flying colours. The spread in the data is limited, and it 
does not noticeably increase with time. The linear kernel (Fig. 
|9j\), on the other hand, shows a different story. Here the mean 
mass, (m), as well as the peak mass, m p - defined as the peak of 
the m 2 f(m) size distribution - evolve exponentially with time. 
Note that the position of m p only depends on the particles that 



contain most of the mass, while (m) is also sensitive to the to- 
tal number of particles. Therefore, (m) lags m p at any time and 
one can show that the gap between the two also increases expo- 
nentially with time. Inevitably, at some point in time, the theo- 
retical value of the m 2 f(m) mass-peak becomes larger than the 
total mass present inside r V. This corresponds to the crossing 
of the dotted line at, e.g., t 20 in Fig.[9j\. In other words, the 
duplication mechanism, needed to conserve N but which has 
the additional effect of enlarging r V, is incapable of keeping up 
with the exponential growth: 'runaway particles' could have 
formed, but the simulated volume *V is just not large enough 
to take them into account. The postulate of the 'duplication 
mechanism' - the particle distribution evolves similarly in and 
outside *V - then breaks down. The only way to avoid this ef- 
fect is to enlarge <V by having more particles in the simulation, 
i.e., to improve the 'numerical resolution'. In Fig.|9j3, we show 
the results, in which A^, instead of N, is held constant (Sect. 
13. 21 . In these simulations N increases with time, starting with 
N = 20000 and ends with more than a million particles. The 
distribution now represents more closely the theoretical curve. 
Growth factors of 14 orders of magnitude in mass (^ 5 in size) 
can then be accurately simulated. 

The drawback, however, is that the computation slows 
down as N increases, since the relative increase in the aver- 
age mass is inversely proportional to N, i.e., A(m)/(m) oc AT -1 . 
These simulations therefore take much more CPU time. It is 
clear that a fundamental limit is reached, in which, given a cer- 
tain CPU power, the calculation of the mass-distributions can 
only be achieved for a limited range. A way to overcome this 
problem is to collide multiple particles per event. In such an 
algorithm collisions are no longer between two particles but 
between two groups of particles. Although this approximates 
the collisional process, the coagulation can be speeded up by 
grouping especially small particles into a single unit. We will 
discuss the grouping algorithm and its implication in a future 
paper. For the simulations in Sect. |4] the product of N and N s 
has been kept constant, which ensures a constant growth (in ex- 
ponential terms) per cycle. We have fixed y/N x N s at 20000 
but made sure that the numerical resolution issues as discussed 
here for the linear kernel, did not occur. Fortunately, realistic 
mass-distributions are not that broad as in the analytic, linear 
kernel; e.g., the smaller particles are quickly removed due to 
Brownian coagulation. 

In summary, we have built an efficient Monte Carlo code 
to follow coagulation. The advantage of this code (above other 
numerical methods) is that it is intuitive, simple to implement 
and expand, and that it takes full account of the stochastic na- 
ture of coagulation. Minor disadvantages are the A^ 2 depen- 
dence on CPU time, the somewhat artificial duplication pro- 
cedure, and the resolution problems shown for the linear ker- 
nel at high m as the simulation progresses. We have addressed 
these concerns and developed methods suitable for the scope 
of this work. DSMC-coagulation methods are very appropri- 
ate to work in conjunction with multi-parameter models. The 
collision model of Sect. [2] with mass, m, and enlargement pa- 
rameter, iff, as variables, can now be put into an evolutionary 
setting. 
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Fig. 9. Tests of the Monte Carlo coagulation code. Linear kernel (Kij = jim, + - (A) fixed N (N = 200 000). The distribution 
function / is computed at times t = (1, 10, 20, 30). At t — 30 insufficient mass is present to provide a good fit to the analytical 
distribution. - (B) fixed N s . To obtain better correspondence with theory, the number of particles, N, is increased as the simulation 
progresses, such that more volume is sampled. N s is fixed at 20000 and N ends with over a million particles. The average 
corresponds well to the theory, yet the amount of CPU time is disproportionally larger at the later times. 



4. Results 

4.1. Application to a protoplanetary disk 

The collision model of Sect.[2]is embedded in the Monte Carlo 
code of the previous section and applied to the protoplane- 
tary disk (Sect. I2.U . The coagulation is restricted to a turbu- 
lent environment in which particles are well mixed with the 
gas, yet do not fragment upon collision. For these reasons, 
we start out with a monodisperse population of sub-micron- 
sized grains (we choose them to be ao = 0.1 /mi) present at 
Z - H g , i.e., at the scaleheight of the disk, where the density is 
a factor exp[-0.5] lower than at the mid-plane. Its collisional 
evolution within the gaseous nebula is followed until the point 
that systematic motions, i.e., settling to the mid-plane, start to 
dominate. Thus, particles are either present and well mixed by 
turbulence, or have started to settle and are no longer in the 
region of interest. Settling is then modelled as a sudden phe- 
nomenon. The reality is, of course, a more gradual transition, 
but the discrete picture here is not a bad ap proximation, since 
the v ertical structure is quickly established dYoudin & Ch iang 
12004 . Settling occurs at the point when the friction time of 
a particle has exceeded a critical friction time, T ra i n , such that 
the scaleheight of the particle, h(Tf), becomes smaller than 
that of the gas, i.e., h(jf) < H g . If self-gravity is neglected 
(valid in the gaseous nebula) and Schmidt numbers, Sc, are 
close to unity, 10 then h can be obtained by equating the par- 



10 The Schmidt number measures the ratio of the gas to parti- 
cle diffusivity; it is supposed to be close to unity if Tf < to. 



ticle diffusion timescale, fdiff = 
tling timescale, f SRt ti = (Q 2 ^ )" 1 



h 2 /vj, with the particle set- 
i.e., Ht£) = H g y/a T /Q.T f (cf. 



umg umescaie, r sett i = \ii Tf ) , i.e., nyrn - n„ yaj/iiTf tci. 
ISchrapler & P enning 2004l lYoudinl20"0 5). The critical friction 



time is then 



M*aiti — n 



(26) 



ISchrapler & Hennin J2004I) 



CgPg <A 2/3 

with a* the compact size of particles (see Fig.|2}. As expected, 
higher values of q-t and higher gas densities (lower Tf) de- 
lay the onset of settling in the sense that the particle has to 
grow further in size before it arrives at the critical friction time. 
Alternatively, increasing if/ also delays settling. 

In the code, settling is implemented as a 'rain-out' : the par- 
ticle is removed from the simulation and the spatial dust density 
decreases. The evolution of these 'rain-out particles' during 
settling is not traced anymore. The focus stays on the particles 
that remain in the layers above the mid-plane. Their evolution 
is followed for 10 7 yr. 

The aj parameter is one of the major uncertainties con- 
cerning the characterisation of protoplanetary disks. One of the 
prime candidates for turbulence is the magneto-rotation al insta- 
bility jHawlev & Balbusll99ltlBalbus & Hawle\ll99ll) . which 
seems to be most robust in well-ionised regions, i.e., in the up- 
per layers of the disk. Another way to chara cterise ar is to re- 
late it to the observed accretion rate, dM/dt, JCuzzi et alJ2 005) 
and then values in the range of 10 4 < aj < 10 2 seem plausi- 
ble. Yet, despite its uncertainty, aj appears in key expressions 
as, e.g., Avij and T ra ; n . Therefore, models are run that cover 
a large range in aj: aj = 10~ 6 - 10~ 2 . Furthermore, we di- 
vide the runs in two categories, which reflect the spatial po- 
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Table 1. Overview of all the runs 



Model-id" 


Model Parameters 


Notes 




R» 


a*r 




6 




RlTa4-P 


1 AU 


10" 


-4 


0.95 


default model 


RlTa2-P 


1 AU 


io- 


-2 


0.95 


increased turbulence 


RlTa6-P 


1 AU 


10" 


-6 


0.95 


decreased turbulence 


RlTa4-C 


1 AU 


io- 


-4 


2/3 


compact model 


R5Ta4-P 


5 AU 


io- 


-4 


0.95 


default model at 5AU 


R5Ta2-P 


5 AU 


io- 


_2 


0.95 


increased turbulence 


R5Ta6-P 


5 AU 


io- 


-4 


0.95 


decreased turbulence 


R5Ta4-C 


5 AU 


io- 


-6 


2/3 


compact model 



" Model names are intended to be mnemonic. 'Rx' stands for ra- 
dius at x AU. Ta, denotes the strength of the turbulence, e.g., 
Ta4 stands for a-f = 10~ 4 . Finally, the suffix denotes whether 
models are porous (P) or compact (C). 

b At 1 AU quartz particles coagulate: p g /pd = 240, y = 
25 ergs cnr 2 , p s = 3.0 g cm 3 ; at 5 AU coagulation is between 
ices: y = 370 ergs cnr 2 , p s = 1.0 g cirr 3 , p g /pd = 57. The gas 
parameters correspond to a minimum mass solar nebula model 
(see Sect.l2~Tl. 

sition in the solar nebula. The 'inner' models correspond to 
conditions at 1 AU where the monomers are quartz with in- 
ternal density p s = 3.0 g cirr 3 and surface energy density of 
y — 25 ergs cirr 2 . The 'outer' models correspond to conditions 
at 5 AU, where the coagulation is that of ices (p s = 1 g cm -3 
and y = 300 ergs cm -2 ) and with an enhanced surfa ce density 
of a f actor 4.2 over the minimum solar nebula ( Naka gawa et alJ 
119861) . For comparison, we also run compact models for the 
aj = 10~ 4 runs. Compact models (denoted by the C-suffix 
in Table are models where the internal structure does not 
evolve, i.e., iff = 1 by definition. An overview of all the models 
is given in Tabled 

4.2. Particle growth and compaction 

In Fig.[TO]mass distributions are shown at various times during 
their collisional evolution. On the y-axis the mass function is 
plotted in terms of m ■ a* ■ f(a*), which shows the mass-density, 
i.e., the mass of grains of compact size a* in logarithmic bins. 
The panels compare the results of compact coagulation (panels 
A and B) with those of porous coagulation (panels C and D) for 
a turbulent strength parameter of aj = 10~ 4 . The coagulation 
is calculated at 1 AU (quartz; panels A and C) and at 5 AU 
(ices, panels B and D). In Fig.[^the a i = 10~ 4 model at 1 AU 
is compared to other a T models at 1 AU (see Table [Q. In Fig. 
[OJaverages of the size distributions of Figs.llOFV, C are shown 
explicitly with time. Here (a) m is the mass-weighted size, 

(a)m = ^r— X (27) 

of the population. Thus, while (a) gives the average particle 
size, (a) m corresponds to the (average) size to which a unit of 
mass is confined. Because of this weighting, (a) m > (a), with 
the equality valid only for monodisperse distributions. 

The porosity evolution is also displayed in Fig. El where 
we plot the ratio of (a) m to (a*) m . (a*) m is defined analogously 



10 




time [yr] 

Fig. 13. Evolution of the size distribution with time. The mass- 
weighted size, (a) m (solid line), and the mean size, (a) (dashed 
line), are calculated for the default models {aj = 10~ 4 ) at 1 AU 
(top lines) and 5 AU (bottom lines). Lines are averaged over 
the 50 simulation runs. After t > 10 3 yr the coagulation drives 
most of the mass into the largest particles. 



to Eq. (I27> . but then for the compact particle size. (a) m /(a*) m 
then gives the mass-weighted averaged enhancement of the ge- 
ometrical radius of the particles (see, Fig. . In the case of 
monodisperse populations or runaway growth one size dom- 
inates the average and the ratio is directly related to (iff) m as 
(a) m /(a) m = (i/f)m This ratio is plotted vs. (a*) m in Figs. [21 
for various models. It shows the build-up of porosity during 
the fractal stage, the stabilisation during the compaction stage, 
upto the stage where the particles rain-out. In Fig. ^1 the av- 
erage sizes of some of these rain-out particles are indicated by 
the detached crosses. Note that (for illustrative reasons) the x- 
axis for these particles corresponds to porous size, a, rather 
than compact size, a*. Properties of the 'rain-out' population 
are given in Table In Fig.[2]the temporal stage is indicated 
by numbered ticks. 

The evolution of the mass-distribution (Figs. llOlll 1111 31 can 
be divided into three stages. Initially, since particles start out 
as grains with sizes of 0.1 fim, Brownian motion dominates. 
The size-distribution is therefore rather narrow, because the 
Brownian collision kernel favours the lighter particles. Quickly 
(~ 10 2 yr), however, turbulent velocities become dominant and 
relative velocities are now highest for the more massive (high 
Tf) particles. Once the first compaction event occurs (dotted, 
blue line 11 ) Tf enters the regime in which it becomes (at least) 
proportional to size and the pace of coagulation strongly ac- 
celerates toward larger sizes. These findings correspond well 

1 1 References to colours only apply to the electronic version of this 
paper. 
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Fig. 10. The mass function plotted at various times for the aj = 10 models. The panels compare the coagulation of the 
compact models (ijf = 1, top panels) with those where porosity effects are included (bottom panels). Left (right) panels show the 
coagulation of quartz (ice) particles at 1 AU (5 AU). Each plot shows the mass function at every logarithmic interval in time from 
t — 10 yr until t = 10 7 yr. In the first ~ 10 2 yr Brownian motion dominates the coagulation. Subsequent evolution is driven by 
turbulence-induced velocity differences and includes the moment of first compaction (blue, dotted line) and first rain-out (red, 
dashed line). After rain-out (f > 10 4 yr), the mass density in the gaseous nebula decreases and the mass function collapses. In the 
compact models the blue, dotted curve also corresponds to the first time that E > E m \\. Grey scales indicate the spread in the 50 
realisations of the simulation. 



with the simple analytical model of Blum| J2004I) . where in 
his Fig. 5 the Brownian motion driven growth is also followed 
by a stage in which the growth is exponential. This evolution 
could be called 'run-away', were it not for the fact that (in our 
case) the particle distribution is truncated at r ra j n (Ea. l26t . The 
distribution at the first rain-out event is shown by the dashed, 
red line. Thereafter, the mass function collapses and evolves 
to a monodisperse population, close to the rain-out size (Figs. 
I10I13I >. Two causes conspire to make these particles favoured: 
first, large particles can only be (efficiently) removed by a colli- 
sion with a similar-sized particle (and no longer by a larger par- 
ticle since these have disappeared); second, in the aj = 10~ 4 
models friction times are always in the tj < t s regime and rel- 
ative velocities between similar-sized particles are suppressed. 



For these reasons, particles in the aj = 10~ 4 models near rain- 
out deplete the smaller particles faster than they deplete them- 
selves, and the size distribution evolves again to monodisper- 
sity. Note, however, that this behaviour is essentially caused by 
the imposed presence of a sharp cut-off size due to the rain-out 
criterion. In reality, a more smooth transition can be expected. 

Although the qualitative trends between the porous and 
compact models are essentially the same - fractal growth, com- 
paction and run-away growth, rain-out and depletion - it is un- 
ambiguously clear that the porous models evolve to larger par- 
ticles, as is also seen in Table |2] in which the values for the 
rain-out particles are tabulated. The size difference at rain-out 
is ~ 2 order of magnitude in size and ~ 4 orders of magnitude 
in mass. Particles only rain-out at T ra i n and in the porous models 



16 C.W. Ormel et al.: Dust coagulation in protoplanetary disks: porosity matters 




10 o 10 10 J l0 z 10 10 u 10 3 10 10 J 10 z 10 10°10 10 J 10 Z 

compact size, a comp [cm] compact size, a comp [cm] compact size, a comp [cm] 



Fig. 11. Effects of turbulence on the coagulation. Panels show the collisional evolution at 1 AU of the porous models, yet with 
ax values of 10~ 2 (A), 10~ 4 (B) and 10~ 6 (C). The scaling of the axis is the same throughout the panels. In the aj = 10~ 2 models 
the spread in the runs is very large, causing the error bars to overlap. In the aj = 10~ 6 model the particles rain-out without 
compacting. 



Evolution curves 1 AU models Evolution curves 5 AU models 
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Fig. 12. Evolution curves of various models. In these panels the enlargement factor is characterised by the ratio of the mass- 
weighted porous size, (a) m , over the mass-weighted compact size, (flcomp)m, ((a*) m in the text), against (a C omp)m- Rising curves 
correspond to an increase of porosity due to fractal growth, and horizontal or declining lines indicate compaction. Numbers give 
the temporal stage of the coagulation (i.e., t = 10')- The detached, coloured crosses indicate the average porous sizes (x-axis) and 
the size enhancement of the rain-out particles (y-axis). 



particles have to grow much further before the critical friction 
time is reached (Ea.l26>. The inclusion of porosity in the coag- 
ulation models thus allows particles (when a T > 10~ 4 ) to grow 
to cm/dm sizes in the gaseous nebula, i.e., before settling to the 
mid-plane. 

Apart from determining the size at which particles rain- 
out, ffj also determines the pace of coagulation. In the aj = 



10~ 2 models (Fig. II coagulation is rapid. Also, a large de- 
gree of stochasticity is seen among different models. In the 
aj = 10~ 6 models, on the other hand, the turbulent veloci- 
ties are very small and the support against gravity is weak, 
such that that rain-out happens before any compaction takes 
place. These models are most reminiscent of the 'laminar 
nebula', where systematic (i.e., settling) velocities dominate 
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model 


{a} 


(m) 


(if/) 




(a/a*) 


Traill 


> 


(r 99% > 




cm 


gr 










yr 




yr 


(1) 


(2) 


(3) 


(4) 




(5) 




(6) 




(7) 


RlTa4-P 


2.1 ± 0.1 


1.4 ± 0.1 


(8.6 


± 0.2) x 10' 


A A 


It u.u 


(4 7 


± 0.2) x 10 3 


(4.5 ± 1.0) x 10 4 


RlTa4-C 


2.4 x 1(T 2 


(1.9 ± 0.2) x 10~ 4 


1 




1 




(2.6 


±0.1)x 10 3 


(4.4 ± 0.3) x 10 4 


R5Ta4-P 


0.72 ± 0.04 


(1.4 ± 0.1) x 10~ 2 


(1.1 


±0.1) x 10 2 


4.8 


±0.1 


(1.9 


±0.1)x 10 4 


(2.2 ± 0.3) x 10 5 


R5Ta4-C 


6.7 x 10~ 3 


1.2 x 10~ 6 


1 




1 




(8.1 


± 0.3) x 10 3 


(2.1 ±0.1) x 10 5 


RlTa2-P 


7.8 ±0.7 


(1.6 ± 0.3) x 10 3 


(3.9 


±1.1) 


1.6 


±0.1 


(2.6 


± 0.3) x 10 2 


(4.0 ± 19) x 10 5 * 


R5Ta2-P 


24 ±2 


(1.7 + 0.3) x 10 3 


(3.6 


±0.1) x 10 1 


3.3 




(1.4 


±0.1)x 10 3 


(4.5 ± 4.6) x 10 3 b 


RlTa6-P 


2.1 x 10~ 2 


1.3 x 10~ 6 


(8.7 


±0.1) x 10 1 


4.4 




(1.3 


±0.1)x 10 3 


2.3 x 10 6 


R5Ta6-P 


1.3 x 10~ 3 


(4.6 ± 0.1) x 10- 10 


(2.0 


±0.0)x 10 1 


2.7 




61 ± 


7.2 


1.1 x 10 5 



8 Entries denote: (1) model-id (see Table0; (2) size at rain-out; (3) mass at rain-out; (4) enlarge- 
ment factor at rain-out; (5) size enhancement at rain-out; (6) time at which first rain-out occurs; 
(7) time-interval over which 99% of the mass has rained-out. Values have been averaged over 
the 50 runs with the error bars reflecting the variations between the 50 runs of the simulation. 
(The spread is only given when the rms-value exceeds the second significant digit.) 



Some simulations did not achieve a 99% rain-out of the density, so that > 10 7 yr. This 
caused the large spread. 



and are therefore most prone to gravitational instability ef- 
fects faubbard & Blackman 2005). The comparison between 
the various ax-models is perhaps best seen in the 'evolution 
tracks' of Figs. ^] They show that initially all porous mod- 
els follow the same (fractal) curve, until the moment com- 
paction occurs. In the 1 AU, aj = 10~ 2 model a significant 
compaction of aggregates can clearly be observed (descending 
line). After t > 10 3 yr, this is followed by a slight increase in 
(a) m /(a comp ) m ; apparently due to the heavy rain-out, most col- 
lisions are again in the fractal regime. 

Comparing the coagulation of the two materials studied 
here, i.e., quartz for the 1 AU models and ice for the 5 AU 
models, one sees similarities in their collisional evolution (see 
also Fig. I13l >. It seems, however, that the coagulation at 5 AU 
is somewhat slower. This can be explained naturally because 
of the lower density, but also perhaps because compaction is 
more difficult to achieve due to the higher surface density (y) 
of ices. Although the differences are subtle, one can see, e.g., 
in Fig.^]that the curves of the aj - 10 -4 , 10~ 2 porous mod- 
els level-off at higher (a) m -to-(a) ratio than the corresponding 
1 AU curves, indicating compaction is achieved 'easier' at 1 
AU. Also, at aj = 10~ 2 the 1 AU particles that rain-out are 
strongly compacted (an even higher aj would have led to frag- 
mentation), while the rain-out particles at 5 AU do not compact 
considerably before rain-out (see also Table|2ji. Thus, the larger 
surface energy (y) of ices translates into a higher rolling energy 
and, subsequently, to decreased compaction. 

5. Discussion 

It now becomes clear that the internal structure of particles, 
represented here by the enlargement parameter iff, is a variable 
of key importance in models of dust-aggregation. To illustrate 
this point further, Fig.^^compares the 'evolution curve' of our 
(default aj = 10~ 4 , 1 AU) model with those of compact, PCA 
and CCA aggregation (grey curves). The superimposed black 
curve connects the (m, ifr) values of the most massive particle 
resulting from the collision model. The small, erratic structure 




10 10 10 5 10° 
mass [g] 

Fig. 14. The m-if/ relation for several aggregation models. 
Plotted is if/(m) for: the most massive particle in one of the 
RlTa4-P models (black curve); the PCA, CCA aggregation 
models, discussed in Sect. 12.41 (grey curves); and compact 
(i/r = 1) models. The dashed lines indicate points of equal fric- 
tion times, T ra i n ~ 500 s and f s w 1 600 s. The cross shows 
the point at which the model experiences the first compaction 
event. 

in this curve corresponds to the fact that the particles fluctu- 
ate stochastically during the simulation. Furthermore, since we 
compare single particles here, instead of a (weighted) mean of 
the distribution as, e.g., in Fig. [2] a fixed point in the figure 
corresponds to one particular friction time and lines of equal 
friction times lie parallel to the dashed lines indicating r = T ra j n 
and t = t s . For the specific choice of aj = 10~ 4 we see that r ra in 
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Fig. 15. Geometrical optical depth to the mid-plane as function of time for the 1 AU models (left) and the 5 AU models (right). 
The optical depth is computed at every logarithmic interval in time and also at the point where particles rain-out (red symbol). 
For illustrative purposes the points are connected by lines. Error bars denote the spread throughout the 50 runs of each model. 



comes prior to f s and velocities therefore remain modest. Still, 
in our collision model the coagulation will reach a point, indi- 
cated by a cross in Fig. ^] at which some collisions, having 
the right mass ratio and relative velocity, are energetic enough 
to cause compaction (E > E ro \\), which effectively halts any 
further increase in porosity. However, due to the earlier exten- 
sive build-up of porosity in the fractal regime, the particle dis- 
tribution now evolves to larger sizes as compared to the com- 
pact models (Fig. 1 1 01 . causing the rain-out masses to be or- 
ders of magnitudes higher (Tabled}- Thus, a major part of the 
growth takes place in the nebula phase. Large, porous particles 
are quickly produced, stay in the nebula mixed with the gas 
and only settle when they are sufficiently compacted, e.g., by 
energetic collisions (this paper) or shocks (see below). 

The curve of our model, with its characteristic bending 
point due to compaction, is a direct result of treating porosity as 
a dynamic variable that is altered by the collisional process. In 
all other models in Fig.^] on the other hand, compaction is not 
incorporated, resulting in straight lines when the growth is in 
the fractal regime. In the PCA/CCA fractal models compaction 
is of course a priori ruled-out. However, the pre-assumed ab- 
sence of compaction in the PCA/CCA models is consistent with 
the low collision energies in these limiting cases. In CCA, fric- 
tion times barely grow (Tf oc m ), the f s 'threshold' is not 
reached and relative velocities vanish for similar particles. In 
PCA, on the other hand, relative velocities are higher, but the 
collision energy is now suppressed by the reduced mass. Thus, 
if the coagulation process would (for some reason) be restricted 
to these limiting growth modes, aggregates will not restructure. 
Fig.ll4lshows how this affects the overall coagulation process. 
It shows, e.g., that compact grains rain-out at ~ 10~ 4 g, 'PCA 
grains' at ~ 10~ 2 g, porosity-evolving particles of our model 
at ~ 1 g, and 'CCA particles' will grow forever! It is clear that 



the porosity evolution of collisional agglomerates is of decisive 
influence on the coagulation process. Modelling the porosity 
evolution in combination with a microphysical collision model 
is therefore a key requirement for a full understanding of the 
first stages of planet formation. 

To quantify the effects of coagulation on the appearance of 
the disk, we have calculated optical depths to the mid-plane in 
the MSN. As a first order approximation, the vertical structure 
is to be taken of constant density and extends over one scale- 
height. We assume that this layer is represented by the parti- 
cles of our simulation and that the rain-out particles (which 
we do not follow) are below it, i.e., at the mid-plane regions 
of the disk. The geometrical optical depth (i.e., at visible/UV- 
wavelengths) to the mid-plane is then calculated as 

1 geom 

= H g J dmf(m)na 2 . (28) 

Results are given in Fig. ^] for the 1 AU and 5 AU models. 
In Fig. I15l\ it is seen that the aj = 10~ 6 models stay optically 
thick for most of the disk's evolution. Note also that the aj — 
10" 4 porous models (solid lines) and the aj = 10 4 compact 
models (dashed-dotted line) do not deviate much in T geom . This 
shows again the dual effects of porosity on the population: it 
increases the geometrical cross-section, yet it also speeds up 
the coagulation, causing more mass to be 'locked' inside large 
particles. At 5 AU (Fig.ll5B). the timescales are longer and the 
disks only becomes optically thin after ~ 10 6 yr when aj < 
10~ 4 . In the aj = 10~ 2 models, evolution to optical thinness is 
very fast at both radii. It is clear that, within the frame-work 
of these models, the inner regions of protoplanetary disks are 
rapidly depleted of small grains unless aj ~ 10~ 6 or less. 

There is a further serious issue hidden here. All our models 
show a rapid decline of the (sub-)micron size population on a 
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timescale of ~ 10 3 yr. This is a well-known problem in mod- 
els for grain coagulation in protoplanetary environments: the 
densities are high e nough for coagulation to proceed rapidly 
JPullemond & Dominikl2005l) . Furthermore, the turbulence in- 
duced relative velocities promote collisions between particles 
with different friction times, i.e., between small and larger par- 
ticles. In contrast, observations reveal the presence of copi- 
ous amounts of small grains in the disk-photospheres of iso- 
lated Herbig AeBe stars and T-Tauri stars, pointing toward 



the presence of small g 


;rains 


on timescales of some 10 yr 


( van Boekel et al. 


2004, 


2005 


: Meeus et al.ll2003l iNatta et alJ 


2006). This discrepancy implies a continuous replenishment of 



(sub)micron-sized grains. In particular, it may reflect the im- 
portance of vaporisation and condensation processes contin- 
uously forming fresh, small grains. Likely, this vaporisation 
and condensation would be localised in the hot and dense, in- 
ner regions of the disk and these grains would then have to 
be transported upwards and outwards to the disk photosphere 
through diffusion processes. The high degree of crystallinity of 
silicates in the inner few AU of protoplanetary disks also points 
toward the i mportance of co ndensation processes in these en- 
vironments (van Boekel et al. 2004). Likewise, the presence of 
crystalline silicates in the cold outer regions of protoplanetary 
disks has been attri buted to large scale mixing of materials i n 
these environments ( Bockelee-Morvan et al. 2002; Gail 2004). 
Alternatively, the replenishment of small grains is through col- 
lisional fragmentation. These energetic collisions could take 
place either in the gaseous nebula due to high relative veloc- 
ities driven by a high aj, or in the mid-plane regions with the 
subsequent upwards diffusion of small grains. While it might 
be difficult to sustain aj > 10 -2 over a prolonged period 
of time, fragmentation in the mid-plane regions seems viable 
since more massive particles will reside here. Furthermore, 
if the mid-plane becomes dust-dominated (pd > p„), shear- 
turbule nce will develop, further augmenting the collisional en- 
ergies JCuzzi et alJl993l) . 

Further constraints on the collisional growth of grains in 
protoplanetary disks are provided by the solar system record; 
specifically, the chondrules and Ca-Al-rich Inclusions (CAI), 
which dominate the composition of primitive meteorites. These 
millimetre-sized igneous spherules are high-temperature com- 
ponents that formed during transient heating events in the early 
solar system. We realise that the cm-sized fluff-balls formed 
in our porous coagulation models are in the right mass-range 
of these meteorite components. It is tempting to identify these 
fluff-balls as the precursors to the chondrules and CAIs. We 
might then envision a model where the collisional evolution 
is terminated by the flash-heating event, for example a shock 
or lightning, which leads to melting and the formation of a 
chondrule and subsequent immediate settling. During the set- 
tling phase the chondrule may acquire a dust rim by sweeping 
up small dust grains or other unprocessed fluff-balls still sus- 
pended in the nebula JCuzzil2004l) . One key point to recognise 
here i s that chondr ules show a spread in age of a few million 
years llWoodl Eb05). which indicates that the collisional grain 
growth process takes place over a much longer timescale than 
our models would predict (see above). 



In this study we have focused on the agglomeration driven 
by random motions - either Brownian or turbulent - high 
up in the nebula. Growth is then presumed to be terminated 
once the aggregate has compacted enough to settle. At that 
point, the aggregate will drop-down about one scaleheight af- 
ter which further growth must occur for further settling to 
continue. In reality, instead of the simple two-component pic- 
ture of nebula and mid-plane presented in this work, the neb - 
ula will acquire a stratified appearance (Dubru lle et alj fl9951. 
where larger particles with higher friction times have smaller 
scaleheights. Collisional evolution models should be able to 
include this stratified nature of the disk. This stratification 
also extends in the radial direction due to radial drift mo- 
tions and turbulent diffusion. Incorporation of these motions 
into the Monte Carlo code will be challenging. We expect 
that the study presented here can serve as the basis for in- 
corporating realistic grain growth in hydrodynamical models, 
likely in the form of well-selected 'collision-recipes'. An obvi- 
ous step would be to include collisions that exceed the Ef mg 
limi t in the collision mod el. Note, for example, that in the 
iDominik & Tielensl i 19971) terminology what we have called 
'fragmentation' should in fact be sub-divided in a continuous 
range of aggregate disruptions. At first monomers will be lost 
and only if E » Ef rag are aggregates completely shattered. 
Other extensions to the model are to allow for a distribution 
of monomer sizes and to use monomers of different chemical 
composition. Both will affect the critical energy for restructur- 
ing, E IO \\. However, with an increasing number of parameters 
characterising the collision, it is worthwhile to verify exper- 
imentally - either through laboratory experiments or through 
detailed numerical calculations - which are of prime impor- 
tance. 

6. Conclusions 

We have presented a model that incorporates the internal struc- 
ture of aggregates as a variable in coagulation models. We used 
the enlargement parameter i// to represent the internal structure. 
It is seen that the internal structure is key to the collisional evo- 
lution, since it crucially affects the dust-gas coupling. However, 
in the model presented in Sect.|2]</' is not a static variable. It is 
altered by the collisional process and we have constructed sim- 
ple recipes to include this aspect in coagulation models. 

Next, we have applied the new collision model to the colli- 
sional evolution of the turbulent protoplanetary disk, until par- 
ticles rain-out to the mid-plane. Our main conclusions are: 

- With the treatment of porosity as a variable, three different 
regimes can be distinguished: fractal growth, compaction 
and fragmentation dBlumll2004 . These regimes are also 
reflected in the collisional evolution of the size distribu- 
tion: fractal growth (mostly Brownian motion), compaction 
(growth accelerates) and rain-out. 

- The collisional evolution of our porosity-evolving model is 
quantitatively different from PCA/CCA aggregation mod- 
els in which porosity can be parameterised by a fixed ex- 
ponent. Therefore, a microphysical collision model is a key 
requirement for coagulation models. 
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- Due to the porous evolution, particles upto dm-sizes can 
be suspended in the gaseous nebula, orders of magnitude 
larger in mass than in models of compact coagulation. 
Therefore, chondrule precursors could have had their ori- 
gin in regions above the mid-plane. 

- If aj < 10~ 2 , no fragmentation occurs in the gaseous neb- 
ula. Therefore, if 10~ 6 < aj < 10~ 2 , the inner nebula will 
become optically thin on timescales of ~ 10 4 yr, unless an 
influx of small grains takes place. 
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Appendix A: List of symbols 



Table A.l. List of symbols. 



Av relative velocity 

E surface density 

Q angular (Keplerian) rotation frequency 

A geometrical cross-section (1 particle) 

C collision rate 

E collision energy 

^restr/^max-c/^frag energy limits for particle restructuring/ 
maximum-compression/ fragmentation 

H g /h gas/particle scale height 

K collision kernel 

N number of monomers (Sect.|2j; number 

of particles (Sect.[3) 

N c number of contacts 

N s number of distinct particles (number of 

species) 

R heliocentric radius 

Re Reynolds number 

T temperature 

V/'V particle/simulated volume 

q-t turbulent viscosity parameter 

£ critical displacement 

y specific surface energy 

S area-mass exponent (A ~ m 5 ) in fractal 
regime 

k number of coagulations 

yu reduced mass 

y ra /v T molecular/turbulent viscosity 

tf/ enlargement parameter 

PilPglPs dust/gas/specific density 

a collision cross-section 

Tf friction time 

Train friction time at which particles are re- 
moved from simulation due to rain-out 

a/ao radius (size)/monomer size 

c g thermal sound speed 

gi number population of species i. 

&b Boltzmann constant 

I eddy length scale 

in mass 

/ particle distribution function (number 

density spectrum) 

r random deviate 

t time 

t s /to turn-over timescale of smallest/largest 
eddy 

v velocity 

v s /v velocity smallest/largest eddies 



